module sw_case8_mod

  use const_mod

  implicit none

  private

  public sw_case8

  real(8), parameter :: radius  = 6371.220d3
  real(8), parameter :: gravity = 9.80616d0
  real(8), parameter :: omega   = 7.292d-5
  real(8), parameter :: umax    = 80.0d0
  real(8), parameter :: lat0    = pi / 7.0d0
  real(8), parameter :: lat1    = pi / 2.0d0 - lat0
  real(8), parameter :: lat2    = pi / 4.0d0
  real(8), parameter :: en      = exp(-4.0d0 / (lat1 - lat0)**2)
  real(8), parameter :: gh0     = gravity * 1.0d4
  real(8), parameter :: ghd     = gravity * 120.0d0
  real(8), parameter :: alpha   = 1.0d0 / 3.0d0
  real(8), parameter :: beta    = 1.0d0 / 15.0d0

contains

  subroutine sw_case8(lon, lat, u, v, gh, ghs, g, f, radi)

    real(8), intent(in ) :: lon
    real(8), intent(in ) :: lat
    real(8), intent(out), optional :: u
    real(8), intent(out), optional :: v
    real(8), intent(out), optional :: gh
    real(8), intent(out), optional :: ghs
    real(8), intent(out), optional :: g
    real(8), intent(out), optional :: f
    real(8), intent(out), optional :: radi

    integer neval, ierr
    real(8) abserr
    
    if(present(radi))radi = radius

    if (present(u  )) u = u_function(lat)
    if (present(v  )) v = 0.0d0
    if (present(gh )) then
      call qags(gh_integrand, -pi0p5, lat, 1.0e-15, 1.0e-10, gh, abserr, neval, ierr)
      gh = gh0 - gh
      gh = gh + ghd * cos(lat) * exp(-(lon / alpha)**2) * exp(-((lat2 - lat) / beta)**2)
    end if
    if (present(ghs)) ghs = 0.0d0
    if (present(g  )) g   = gravity
    if (present(f  )) f   = 2 * omega * sin(lat)

  end subroutine sw_case8

  real(8) function u_function(lat) result(res)

    real(8), intent(in) :: lat

    if (lat <= lat0 .or. lat >= lat1) then
      res = 0
    else
      res = umax / en * exp(1.0d0 / (lat - lat0) / (lat - lat1))
    end if

  end function u_function

  real(8) function gh_integrand(lat) result(res)

    real(8), intent(in) :: lat

    real(8) u, f

    u = u_function(lat)
    f = 2 * omega * sin(lat)
    res = radius * u * (f + tan(lat) / radius * u)

  end function gh_integrand

end module sw_case8_mod
